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ABSTRACT 

, The dynamics of two massive black holes in a rotationally supported nuclear disc 

of mass Mdi S c = 1O 8 M is explored using N-Body/SPH simulations. Gas and star 
particles are co-present in the disc. Described by a Mestel profile, the disc has a 
£N| ■ vertical support provided by turbulence of the gas, and by stellar velocity dispersion. 

t> ' A primary black hole of mass 4 x 10 6 M Q is placed at the centre of the disc, while a 

, secondary black hole is set initially on an eccentric co-rotating orbit in the disc plane. 

Its mass is in a 1 to 1, 1 to 4, and 1 to 10 ratio, relative to the primary. With this 
choice, we mimic the dynamics of black hole pairs released in the nuclear region at 
the end of a gas-rich galaxy merger. It is found that, under the action of dynamical 
friction, the two black holes form a close binary in ~ 10 Myrs. The inspiral process is 
\ insensitive to the mass fraction in stars and gas present in the disc and is accompanied 

-»■»« . by the circularization of the orbit. We detail the gaseous mass profile bound to each 

' black hole that can lead to the formation of two small Keplerian discs, weighing « 2% 

of the black hole mass, and of size ~ 0.01 pc. The mass of the tightly (loosely) bound 
O ' particles increases (decreases) with time as the black holes spiral into closer and closer 

\ orbits. Double AGN activity is expected to occur on an estimated timescale of < 10 

c/3 , Myrs, comparable to the inspiral timescale. The double nuclear point like sources that 

may appear during dynamical evolution will have typical separations of < 10 pc. 

Key words: black hole physics - hydrodynamics - galaxies: starburst - galaxies: 
. evolution - galaxies: nuclei 
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1 INTRODUCTION 

Collisions of spiral galaxies may lead to the formation of 
(ultra-)luminous infrared galaxies (U-LIRGs). A clear ex- 
ample is NGC 6240, a ULIRG composed of two gas-rich 
galaxies on the verge of merging, hosting two spatially dis- 
tinct X-ray sources. Buried in a wide-spread star burst, 
these two hard X-ray sources have been identified as a pair 
of accreting massive black holes (MBHs) (Komossa et al. 
2003; Risaliti et al. 2006). 

The majority of (U-)LIRGs are unsettled remnants of 
colliding galaxies undergoing episodes of intense star for- 
mation (see for a review Sanders & Mirabel 1996). A large 
number of (U-)LIRGs hosts a central rotationally supported 
massive gaseous disc (up to 1O 1O M0) extending on scales 
of ~ 100 pc (Sanders & Mirabel 1996; Scoville, Yun & 
Bryant 1997; Downes & Solomon 1998; Tacconi et al. 1999; 
Bryant & Scoville 1999; Greve et al. 2006). As shown in 



numerical simulations, these discs may be the end-product 
of gas-dynamical and gravitational torques (excited during 
the merger), driving large amounts of gas in the core of the 
remnant, where the MBHs are expected to reside (Barnes & 
Hernquist 1991, 1996; Kazantzidis et al. 2005). Inside these 
massive self-gravitating discs, MBH pairs continue their dy- 
namical evolution, and can accrete gas. 

The formation and evolution of MBH pairs in merg- 
ing galaxies have a large number of potentially interesting 
astrophysical consequences. In gas-rich mergers, MBH bi- 
naries can produce, beside double active nuclei, other pe- 
culiar features such as periodic modulations (as seen in OJ 
287), and wiggling jets (see Komossa 2006 for a review). In 
dry mergers, binary MBHs can erode the inner stellar cusp 
(Milosavljevic et al. 2002), creating a population of hyper- 
velocity stars in galactic halos (Yu & Tremaine 2003; Brown 
et al. 2006). Ultimately, binary MBHs may eventually coa- 
lesce emitting low-frequency gravitational waves detectable 
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by LISA (Bender et al. 1994; Hughes 2002; Sesana et al. 
2005). 

Mayer et al. (2006) have shown that in dissipative merg- 
ers, large-scale inflows ending in the formation of massive 
turbulent nuclear discs facilitate the pairing of the MBHs 
hosted in the parent galaxies, and that an eccentric Keple- 
rian binary forms. In this paper, we will explore the phase 
following the formation of the binary. In a previous work 
(Dotti, Colpi & Haardt 2006, hereafter DCH06), we stud- 
ied the dynamical evolution of MBH pairs orbiting inside a 
massive gaseous disc (see also Escala et al. 2005), exploring 
different MBH mass ratios and initial orbital eccentricities. 
We showed that, during MBH inspiral, gas-dynamical fric- 
tional forces circularize the orbits (if initially corotating with 
the disc), and only after circularization a sizable amount of 
gas is captured by each MBH. The force resolution in such 
early simulations was not sufficient to resolve the small ac- 
cretion discs that form around the inspiralling MBHs. 

In the simulations presented in DCH06, the gaseous disc 
was embedded into a larger scale, spherically symmetric stel- 
lar distribution. Fragmentation of the disc was prevented by 
using an adiabatic equation of state, and star formation was, 
simply, neglected. Thus, these simulations did not consider 
the possible co-presence of an axisymmetric nuclear stellar 
disc. Indeed, nuclear stellar discs are observed in many early 
and late type galaxies (van den Bosch et al. 1998; Scorza & 
van den Bosch 1998; Morelli et al. 2004; Krajnovic & Jaffe 
2004; Lopes et al. 2006; Mueller Sanchez et al. 2006; Davies 
et al. 2007). Nuclear stellar discs result from galaxy merg- 
ers and/or from the secular evolution of bars (Krajnovic & 
Jaffe 2004; Kormendy et al. 2005; Ferrarese et al. 2006). 
Davies et al. (2007), observing young stellar discs in NGC 
1068 and NGC 1097, estimate masses of ~ 10 8 M Q , radii of 
~ 50 pc, and a vertical scale heights of ~ 5 — 10 pc, sug- 
gesting that the star velocity dispersion provides significant 
pressure support. 

Aiming at a better understanding of i) the role of a nu- 
clear stellar disc in driving the orbital evolution of the MBH 
binary, and ii) the formation of an accretion disc around each 
MBH, we run a new series of simulations employing the par- 
ticle splitting technique (Kitsionas & Whitworth 2002), al- 
lowing us to improve the numerical resolution with respect 
to DCH06 and other similar simulations (Kazantzidis et al. 
2005). In the most accurate simulation, we can resolve, on 
sub-parsec scales, the gravitational sphere of influence of 
each MBH. For the first time, we detail the mass profile of 
the gas bound to the MBHs, so that we can conjecture on 
the accretion process during the MBH pairing and inspiral. 

The paper is organized as follows. In Section 2 we de- 
scribe the equilibrium structure of the gaseous and stellar 
disc, and the initial conditions for our different runs. In Sec- 
tion 3 we present the results of our simulations carried on 
varying the MBH masses and the stellar content in the disc, 
as well as the aforementioned simulation at higher resolu- 
tion, and finally, in Section 4 we draw our conclusions. 



2 SIMULATION SETUP 

We follow the dynamics of MBH pairs in nuclear discs us- 
ing numerical simulations run with the N-Body/SPH code 
GADGET (Springel, Yoshida & White 2001). 



In our models, a primary MBH of mass Mi is placed at 
the centre of a gaseous and/or stellar disc, embedded in a 
larger scale stellar spheroid. A secondary MBH of mass M2 
is placed on a bound orbit in the disc plane. 

The gaseous disc is modeled with 235331 particles, has 
a total mass Mnisc = 10 8 Mq, and follows a Mestel surface 
density profile 

S(i?)=S (^), (1) 

where R is the radial distance projected into the disc plane, 
and Eo is the density at scale radius Ro- Fig. 1 shows the 
profile E(-R), in physical units. The disc is rotationally sup- 
ported in R, has a finite radial extension of 100 pc, and 
finite vertical thickness of 10 pc. Initially, the gaseous par- 
ticles are distributed uniformly along the vertical axis. The 
SPH particles evolve adiabatically, and the initial internal 
energy density profile scales as: 



u(R) = KR 



-2/3 



(2) 



where K is a constant defined so that the Toomre parameter 
of the disc, 



Q = 



ttGE' 



(3) 



is > 3 everywhere, preventing disc fragmentation and for- 
mation of large scale over-densities, such as bars and spiral 
arms. In equation (3) k is the local epicyclic frequency, and 
c s the local sound speed of the gas, that, for our choice of 
K, is w 30kms _1 at R = 50 pc. The internal energy of gas 
particles mimics the internal, unresolved turbulence, and, as 
a consequence, Ca has to be considered as the local turbulent 
velocity. 

The spheroidal component (bulge) is modeled with 10 s 
collisionless particles, initially distributed as a Plummer 
sphere with mass density profile (shown in Fig. 2) 



, , 3 MBuigc f, . r 2 



-5/2 



(4) 



where b (= 50 pc ) is the core radius, r the radial coordinate, 
and MBui g c(= 6.98MDi sc ) the total mass of the spheroid. 
With such choice, the mass of the bulge within 100 pc is 
5 times the mass of the disc, as suggested by Downes & 
Solomon (1998). 

The primary MBH has a mass Mi = 4 x 1O 6 M , while 
for the secondary we consider three different values, equal 
to 4 x 1O 5 M , 1O 6 M and 4 x 1O 6 M , respectively. The 
secondary MBH is placed inside the massive disc on an ec- 
centric orbit with eccentricity e ~ 0.7, at a separation of 50 
pc from the central MBH. 

We evolve our initial composite model (bulge, disc and 
primary MBH) for « 3 Myrs, until the bulge and the disc 
reach equilibrium. Given the initial homogeneous vertical 
structure of the disc, the gas initially collapses on the disc 
plane exciting small waves that propagate through the sys- 
tem. The spheroid is gravitationally stable on large scales 
(r > 20 pc) and contracts in the central region until equilib- 
rium is reached. Fig. 1 and 2 show the equilibrium profiles, 
and in particular the density enhancement in the central re- 
gion of the spheroid. The stellar density increase does not 
affect the inner region (r < 10 pc) of the (denser) disc nor 
the dynamics of the MBHs. 



© 0000 RAS, MNRAS 000, 000-000 



Binary black holes 3 



10 6 p 




R [pc] 

Figure 1. Gaseous disc surface density £ as a function of distance 
R from the central MBH. Red dashed line refers to the surface 
density of the Mestel disc (i.e., the disc before it reaches equi- 
librium), while black solid line describes the equilibrium profile 
after an elapsed time rj 3 Myrs. 



10 4 p 




r [pc] 

Figure 2. Density p of the stellar spheroid as a function of dis- 
tance r form the central MBH. Red dashed line refers to the initial 
Plummer profile, while black solid line refers to the equilibrium 
profile. 



We run four different sets of simulations (for a total of 
12 according to the mass M2 of the secondary MBH), assum- 
ing a purely gaseous disc, and a disc in which 1/3, 2/3, and, 
finally, all gas particles are turned into collisionless particles, 
respectively. For each disc model with fixed star fraction, we 
evolved the initial condition in isolation until equilibrium is 
reached, as indicated above. The gas — > star conversion in 
the disc is obtained converting randomly selected gas parti- 
cles in collisionless particles. The latter own the same mass 
and position of the initial gas particles, so that the stellar 
disc component follows the same density distribution of the 
gas. The velocities of the transformed particles are modi- 
fied by adding an isotropic component equal (in modulus) 
to the local sound speed, preventing the vertical collapse of 
the new disc. Gas particles are transformed at a constant 
rate of « 16OM0yr -1 . This rate prevents spurious relax- 
ation that could, in principle, change the structure of the 
disc. The outputs of this trial simulation at three different 
times (t « 0.2, 0.4, and 0.6 Myrs) correspond to a frac- 
tion of 1/3, 2/3, and to the total disc mass transformed in 
stars. These outputs are used as new initial conditions for 
the three sets of simulations with an axisymmetric stellar 
population. We do not convert any gaseous particle in stars 
when we follow the dynamics of the MBHs, so that the disc 
stellar fraction remains constant. The secondary MBH is in- 
serted in the plane of the disc, after the composite system 
has relaxed to equilibrium. 

The number of particles used to model the two compo- 
nents allows for a spatial resolution ~ 1 pc, corresponding 
to the gravitational softening which is the same for the colli- 
sionless, gaseous and MBH particles. The thermodynamical 
parameters of the gaseous disc are evaluated averaging over 
a subsample of 50 neighbours. 

The mass of the primary and the mass ratio q = 
M2/M1 fix two characteristic scale-lengths. Following Mer- 
ritt (2006), we define the gravitational influence radius of 
the larger MBH (rinf) as the radius of a sphere that encloses 
a mass (in gas and stars) equal to 2Mi. For separations 
smaller than r- m f the two MBHs form a "binary". For the 
Plummer sphere and Mestel disc parameters used, n n f w 6 
pc, so it is well resolved in our simulations. When the binary 
forms, we can estimate the hardening radius ah. Assuming 
a SIS profile, a h = 0.25gr inf /(l - g) 2 (Merritt 2006), and, 
for q = 1, a h » 1 pc, while for q = 0.25 (0.1) a h « 0.3 (0.1) 
pc. 

The main input parameters of our simulations are sum- 
marized in Table 1. We also performed a simulation at higher 
resolution, starting from an intermediate snapshot of run 
Al, that will be detailed and discussed in section 3.2. 



3 RESULTS 

3.1 Orbital decay and circularization 

Fig. 3 shows the MBH relative separation s and eccentricity 
as a function of time (upper and lower panel respectively) 
for equal mass binaries (runs Al, Bl, CI, and Dl). Regard- 
less the fraction of star-to-gas disc particles, the secondary 
MBH spirals in at the same pace, reaching the force resolu- 
tion limit of ~ 1 pc after ~ 5 Myrs. The velocity dispersion 
of stars is similar to the gas sound speed, and the two com- 
ponents share the same differential rotation. This is why the 



© 0000 RAS, MNRAS 000, 000-000 



4 Dotti, Colpi, Haardt, Mayer 



Table 1. Run parameters 



run 








Mi. 

Disc 


Mi , 

Bulg 


e 

e 


Al c 






4 








A 2 





1 


1 


100 


698 


0.7 


A3 






0.4 








Bl 






4 








132 


1/3 


4 


1 


100 


698 


0.7 


B3 






0.4 








CI 






4 








02 


2/3 


i 


1 


100 


698 


0.7 


C3 






0.4 








Dl 






4 








D2 


1 


4 


1 


100 


698 


0.7 


D3 






0.4 









a /* : disc mass fraction in stars. 
b Masses are in units of 10 6 Mq. 

c Simulation Al was run also using the particle splitting tech- 
nique to improve force resolution (down to ~ 0.1 pc, see sect. 
3.2). 



i | , , , | , , , r 




t [10 6 yr] 

Figure 3. Equal mass MBHs. Upper panel: separations s (pc) 
between the MBHs as a function of time. Lower panel: eccentricity 
of the MBH binary as a function of time. Black, blue, red and 
green lines refer to stellar to total disc mass ratio of 0, 1/3, 2/3 
and 1 (run Al, Bl, CI, and Dl) respectively. 



dynamical friction on M2 caused by stars and gas is similar. 
As the orbit decays, the eccentricity decreases to e < 0.2. 
This value is not a physical lower limit, but rather a numer- 
ical artifact due to the finite resolution. Circularization was 
found in DCH06 for a full gaseous disc. Here, we show that 
circularization occurs regardless the nature of the disc parti- 
cles (gas and/or stars). Note that circularization takes place 
well before the secondary feels the gravitational potential of 
Mi, so the MBH mass ratio does not play any role in the 
process. 

To show how the circularization process works, consider 
run CI, where both stellar and gaseous particles are present 
in the disc. In Fig. 4 we plot the gas (left panels) and stellar 
(right panels) densities in the disc at two different times, cor- 
responding to the first passage at pericentre (upper panels), 
and at apocentre (lower panels). The green curve shows the 
counterclockwise corotating orbit of the secondary MBH. 
Near to pericentre, the MBH has a velocity larger than the 
local rotational velocity, so that dynamical friction causes 
a reduction of the velocity of the MBH. A wake of parti- 
cles lags behind the MBH trail. On the other hand, near 
to apocentre, the MBH velocity (mainly tangential) is lower 
than the disc rotational velocity, and, in this case, the force 
increases the MBH angular momentum: the wake is dragged 
in front of the MBH trail. Thus, as a result of differential 
rotation, the wake reverses its direction at apocentre accel- 
erating tangentially the MBH, leading to circularization of 
M2 orbit (see also Fig. 3). This seems a generic feature, 
regardless the disc composition, as long as the rotational 
velocity exceeds the gas sound speed and the stellar veloc- 
ity dispersion, as in all cases explored. We remark that in 
spherical backgrounds, dynamical friction tends to increase 
the eccentricity, both in collisionless (Colpi, Mayer & Gover- 
nato 1999; van den Bosch et al. 1999; Arena & Bertin 2007) 
and in gaseous (Sanchez-Salcedo & Brandenburg 2001) en- 
vironments. 

Fig. 5 shows the MBH separation and eccentricity evo- 
lution for the cases with M2 = 10 6 Mq (runs "2" for all sets). 
The color coding is the same as in Fig. 3. As the dynami- 
cal friction time-scale is tx 1 /M2 (for given disc properties) , 
in this case M2 orbital decay occurs on a time » 4 longer 
than in runs "1". We can also observe that, in the case of 
a pure gaseous disc (run A2), circularization is more effi- 
cient compared to all other runs. When the secondary MBH 
moves along the eccentric orbit, its motion is supersonic for 
nearly all the entire orbit. Since for a supersonic perturber 
in a gaseous background dynamical friction is ~ 2 times 
more efficient that in a stellar environment (Ostriker 1999), 
the circularization efficiency is increased (while, as discussed 
above, the orbital decay time-scale is nearly independent on 
the nature of the disc particles) . 

The same scalings hold for runs "3", where M% = 
4 x 10 5 Mq . The differences in the orbital decay timescale 
among runs A3, B3, C3 and D3, are however larger than 
what found in runs "1" and "2", because of the smaller 
secondary-to-background particles mass ratio. In runs "3" 
a random component of the motion of M2 exists, mainly 
caused by the stellar bulge. In fact, the mass ratio between 
M2 and a single bulge particle is 57, so that the background 
can not be treated as a smooth fluid. 
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Figure 4. Snapshots from run CI. The panels show a face-on pro- 
jection of the disc and MBH positions at two different times. The 
color coding indicates the z-averaged gas density (left panel) and 
star density (right panel) in logarithmic scale (units: Mq/pc 3 ). 
The green line traces the MBH counterclockwise prograde orbit. 
The gaseous (stellar) disc radius is R 100 pc. 




2 4 6 8 10 12 14 16 
t [10 6 yr] 



Figure 5. Same as Fig. 3, but for MBH mass ratio 1/4 (runs "2" 
in all sets). Upper panel: separations s (pc) between the MBHs as 
a function of time. Lower panel: eccentricity as a function of time. 
Black, blue, red and green lines refer to stellar to total disc mass 
ratio of 0, 1/3, 2/3 and 1 (run A2, B2, C2, and D2), respectively. 



3.2 Nuclear discs around the spiraling MBHs: the 
high resolution run 

We run a higher resolution simulation to study the eccen- 
tricity and orbital evolution on scales lower than 1 pc. The 
new initial condition is obtained resampling the output of 
run Al (for equal mass MBHs) with the technique of particle 
splitting. Resampling is performed when the MBH separa- 
tion is ~ 14 pc (corresponding to ~ 4 Myrs after the start 
of the simulation). We split each gas particle into AT c h — 10 
"child" particles (Kitsionas & Whitworth 2002), randomly 
distributed around the position of the original parent par- 
ticle within a volume of size ~ hp, where h p is the gravita- 
tional softening (of the parent particles). The velocities and 
temperatures of the child particles are set equal to that of the 
parent ones. Each child particle has a mass equal to 1/N c h 
of the mass of the parent. Splitting is applied to all particles 
whose distance from the binary centre of mass is < 42 pc, so 
that the total number of particles increases only by a factor 
~ 4, while the local mass resolution in the split region is 
comparable to that of a standard ~ 2 x 10 6 particle simu- 
lation with uniform resolution. Our choice of the maximum 
distance for splitting is conservative, aimed at preventing 
that more massive, unsplit gas particles reach the binary 
on a timescale shorter than the entire simulation time. In 
the central split region, the high mass resolution achieved 
fulfills the Bate & Burkert (1997) criterion for gravitational 
softening values down to 0.1 pc. In other words, the resolu- 
tion scales of hydrodynamical forces (SPH smoothing) and 
gravitational forces (gravitational softening) are similar. In 
particular, a sphere with radius « 0.1 pc (corresponding to 
the new gravitational softening) centred on a particle con- 
tains > 1 SPH kernel (.A ne igh = 50 particles for our simula- 
tions) . In Fig. 6 we compare the surface density profile of the 
circum-nuclear gaseous disc in run Al at t = 4 Myrs, for the 
low and high resolution cases. The two profiles differ only 
below the scale of the low resolution limit R < 3 pc. The 
decrease of the gravitational softening corresponds to intro- 
ducing a deeper potential well of the MBH within a sphere 
defined by the former softening radius. Therefore, with the 
improved resolution, the central surface density increases as 
the gas reaches a new hydrostatic equilibrium closer to the 
MBH, as shown in Fig. 6 (red line). The lack of noticeable 
differences in the surface profile at separations R > 3 pc 
confirms the accuracy of the particle splitting technique. 

Results of the high resolution run are shown in Fig. 7. 
The separation decays down to 0.1 pc in ~ 10 Myrs. We 
notice that the ellipsoidal torque regime, defined in Escala 
et al. 2004 (see also Escala et al. 2005) and characterized 
by a fast orbital decay, is not present. We attribute this dif- 
ference to the different thermal state of the gas surrounding 
the MBHs. For a comparison, we rescale the results of our 
simulation to the parameters used in Escala et al. (2005). 
After rescaling, our initial disc is found to be hotter by a 
factor » 2.4 than the hottest disc model presented in Es- 
cala et al. (2005). Since the efficiency of angular momentum 
transport by ellipsoidal deformations decreases with increas- 
ing temperature, in our hotter disc dynamical friction driven 
torques still dominate over the ellipsoidal torques. 

In the high resolution run, the dynamical evolution of 
the MBHs is initially identical to the low resolution case, 
as shown in Fig. 7. Because of particle splitting, the system 
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R [pc] 

Figure 6. Surface density profile of the circum— nuclear gaseous 
disc in run Al at t = 4 Myrs. Black line refers to the surface 
density in the low resolution simulation, red line refers to the 
high resolution (split) simulation. The dashed vertical line marks 
the resolution limit in the un— split simulation. 
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Figure 8. Mass of bound gas particles as a function of binary 
separation s in the final stages (from sa 6 to fa 8.5 Myrs) of the 
high resolution simulation (split run Al). Black lines refer to gas 
particles bound to the MBH initially at rest (primary), red lines 
refer to gas particles bound to the secondary MBH. Solid, dashed, 
and dot— dashed lines refer to WB, B, and SB particles. 
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Figure 7. Separations s (pc) between the MBHs as a function of 
time. Red (black) line refers to the (un-)split run Al. The dashed 
area highlight when the MBHs separation is < 1 pc (correspond- 
ing to the low numerical resolution) in the split run Al. 



granularity is reduced, and therefore the force resolution in- 
creases. In the high resolution run, the binary decreases its 
eccentricity to ~ (before the new spatial resolution limit 
is reached). 

A resolution of 0.1 pc allows us to study the properties 
of the gas bound to each MBH. To this purpose, it is useful to 
divide gaseous particles, bound to each MBH, into three sub- 
sets, according to their total energy relative to the MBH. We 
then define weakly bound (WB), bound (B), and strongly 
bound (SB) particles according to the following rule: 

C (WB) 
E < I 0.25 W (B) (5) 
[ 0.51V (SB), 

where E is the sum of the kinetic, internal and gravita- 
tional energy (per unit mass), the latter referred to the 
gravitational potential W of each individual MBH. Here- 
after WBPs, BPs, and SBPs will denote particles satisfying 
the WB, B, or SB condition, respectively. Note that, with 
the above definition, SBPs are a subset of BPs, which in 
turn are a subset of WBPs. 

We find that the mass collected by each MBH, rela- 
tive to WB, B, and SB particles is JV/wbp ~ 0.85Mmbh ~ 
3.4 x 1O 6 M , Af B p ~ 0.41M M bh ~ 1.6 x 10 6 M Q , and 
M S bp ~ 0.02M M bh ~ 8 x 10 4 M Q , respectively (here 
Mmbh = Mi — M2 of run Al). These masses are of the 
same order for both the primary and secondary MBH, and 
remain constant as long as the MBH separation is s > 1 
pc. At shorter separations WBPs and BPs are perturbed by 
the tidal field of each MBH, and at the end of the simula- 
tion Mwbp and Mbp are reduced as shown in Fig. 8. During 
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Figure 9. High-resolution simulation at t = 4.5 Myrs. Upper 
panel: gaseous cumulative mass fraction of bound particles (for 
the SBPs, BPs, and WBPs, respectively) as a function of the dis- 
tance from each MBH. Lower panel: z— component of cumulative 
orbital angular momentum, per unit bound mass, as a function of 
the distance from each MBH. Black (red) lines refer to the initial 
central primary (orbiting secondary) MBH. Solid, dashed, and 
dotted— dashed lines refer to WBPs, BPs, and SBPs, respectively. 



the same period of time, Msbp associated to the primary 
(secondary) MBH increases by a factor » 4 (w 2.5). This 
result is unaffected by numerical noise since the number of 
bound particles (associated to each class) is > 1 SPH kernel 

(iVnoigh = 50). 

The radial density profiles of WBPs, BPs, and SBPs 
are well resolved during the simulation. The upper panel of 
Fig. 9 shows the normalized mass within distance r from 
each MBH, for WBPs, BPs, and SBPs, at t = 4.5 Myrs, 
when the binary separation is ~ 5 pc. The lower panel 
of Fig. 9 shows the specific angular momentum perpen- 
dicular to the disc plane (L z ) of the particles, computed 
in each MBH frame. Bound particles have a net angular 
momentum with respect to each MBH, and form a pres- 
sure supported ellipsoid. The half-mass radius is similar 
for the two MBHs: ~ 3 pc, ~ 1 pc, and ~ 0.2 pc for 
WBPs, BPs, and SBPs, respectively. The disc gas density 
can be as high as 10 r cm~ 3 . It is then conceivable that, 
at these high densities, dissipative processes could be im- 
portant, possibly reducing the gas internal (turbulent and 
thermal) energy well below the values adopted in our simu- 
lations. If any dissipative precess could reduce efficiently the 
gaseous internal energy, we expect that the bound gas will 
form a cool disc with Keplerian angular momentum compa- 
rable to what we found in our split simulation (see Fig. 9). 
Since L z = U G Mmbh i?MBH,disc, we obtain, for the primary 
MBH, an effective radius 7?MBH,disc ~ 0.1 pc (0.03 pc) for 
WBPs (BPs). The secondary MBH is surrounded by parti- 
cles with a comparatively higher angular momentum (see red 



lines in Fig. 9, lower panel), with a corresponding effective 
radius i?MBH,di ac « 1 (0.13) pc for WBPs (BPs). Finally, 
for SBPs, both holes have 7?MBH,disc <S 0.01 pc, which is 
more than an order of magnitude below our best resolution 
limit. These simple considerations indicate that a more real- 
istic treatment of gas thermodynamics is necessary to study 
the details of gas accretion onto the two MBHs during the 
formation of the binary, and the subsequent orbital decay. 
Nonetheless, our simplified treatment allows us to estimate 
a lower limit to the accretion timescale, assuming Eddington 
limited accretion: 

tacc = r Ed d In ( 1 + ,^ acc ) , (6) 

l — £ \ MmbH.O / 

where e is the radiative efficiency, TEdd is the Salpeter time, 
Mmbh.o is the initial MBH mass, and M acc is the accreted 
mass. Assuming e = 0.1, Eddington limited accretion can 
last for ~ 30 Myrs, ~ 15 Myrs, and < 1 Myrs, if the MBHs 
accrete all the WBPs, BPs, and SBPs, respectively. 



4 DISCUSSION 

In this paper we have shown that dynamical friction against 
a gaseous and/or stellar background is responsible for the 
inspiral of MBH pairs inside massive nuclear discs. A MBH 
binary forms at a separation ~ 5 pc. The density and ther- 
mal distributions of the gas and star particles around the 
binary MBH are such that dynamical friction keeps acting 
on each MBH down to a separation of ~ 1 pc since the co- 
herence of the density wakes excited by the MBHs is main- 
tained down to this scale. Thanks to the particle splitting 
technique applied to the case of a purely gaseous disc, we fol- 
low the binary orbital decay down to « 0.1 pc, and we show 
that angular momentum losses by friction, along co-rotating 
orbits, reduce the orbital eccentricity to a value consistent 
with zero. 

Other processes neglected in our study can help to 
shrink the binary orbit. For example, three-body encounters 
with background stars may become important at MBHs sep- 
arations < rjnf ~ 6 pc (Quinlan 1996; Milosavljevic & Mer- 
ritt 2001; Merritt 2006). The cumulative effect of this colli- 
sional process, studied mainly in spherical backgrounds, can 
lead to an increase of the eccentricity, thus acting against 
the circularization driven by the large-scale action of the 
gaseous and/or stellar disc (Quinlan 1996; Aarseth 2003; 
Matsubayashi, Makino & Ebisuzaki 2005; Berczik et al. 
2006; Sesana et al. in preparation). The effect of three body 
encounters on the MBH orbits when stars form a rotation- 
ally supported disc is still unclear. 

We also quantified the structure of gas particles which 
become bound to each MBH during the orbital evolution. 
We found that during the orbital inspiral, a gas mass ~ 50% 
of the MBH mass is conveyed inside the MBH sphere of in- 
fluence, and that a gas mass m 2% of the MBH mass binds 
deeply to each single MBH. The radial distribution of the 
most bound gas particles suggests that an active Eddington- 
limited accretion phase may set in, for a time < 1 Myrs 
around both MBHs. Since we have neglected gas cooling in 
the simulation, this mass likely represents a lower limit. The 
active phase could last for a longer time (> 10 Myrs), com- 
parable to the inspiral timescale, if all the bound mass is 
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accreted. This highlights the possibility of revealing double 
AGN activity, on spatial scales s < 10 pc. However, since we 
expect that star formation is still ongoing in the disc while 
the MBHs spiral in, double AGN activity could be heavily 
obscured. In the same high resolution simulation, we mea- 
sured the cumulative angular momentum (perpendicular to 
the disc) of particles bound to each individual MBH, al- 
lowing an estimate of the size of the Keplerian disc that 
is expected to form. We found disc sizes of ^ 0.1 pc and 
<C 0.01 pc for the bound and the strongly bound gas com- 
ponents, respectively. Such small extension of the accretion 
disc around each MBH could preserve the gas against tidal 
perturbation and stripping, at least until the binary reaches 
separations of the same order. 

We plan to carry on higher resolution simulations with 
more realistic input physics and radiative cooling, in order to 
follow the last phases of the gas-dynamical evolution of the 
nuclear discs. In particular, we aim at tracing the expected 
transition of the two discs into a single, circumbinary disc 
surrounding the two MBHs. 

Armitage & Natarajan (2002) have shown that the in- 
teraction between the two MBHs and the circumbinary disc 
can drive the binary to coalescence in ~ 10 Myrs, for MBH 
separations < 0.1 pc. They have also shown that, during the 
process, the MBH binary orbital eccentricity increases (Ar- 
mitage & Natarajan 2005). Our planned future simulations 
will test in a self-consistent way such predictions, and will 
explore the possibility of exciting flaring activity during the 
latest phases of MBH binary inspiral. This will constrain the 
properties of an electromagnetic precursor associated to the 
binary coalescence (Armitage & Natarajan 2002; Kocsis et 
al. 2005; Milosavljevic & Phinney 2005; Dotti et al. 2006). 
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